% This scripts create synthetic EEJ time series based on the time stamps
% and longitudes of actual CHAMP equatorial crossings, using know amplitude
% and phases for each of the major tidal components and invert them using
% different tidal modes. The idea is to study the aliasing of spectrum from
% one frequency/ wavenumber to another
%% Load data

load  /Users/manojnair/projects/tides/eef_data_2000_2010.mat;

data_max = length(eej_data_2000_2010);

% All data
% date_min = min(eej_data_2000_2010(:,1));
% date_max = max(eej_data_2000_2010(:,1));

% Aug 2000 - Aug 2005
% date_min = min(eej_data_2000_2010(:,1));
% date_max = datenum(2005,08,1);
% 
% Aug 2005 - Aug 2010
date_max = max(eej_data_2000_2010(:,1));
date_min = datenum(2005,08,1);

L = eej_data_2000_2010(:,1) >= date_min & eej_data_2000_2010(:,1) <= date_max;

eej_data_2000_2010(L,:) = [];

datemat = datevec(eej_data_2000_2010(:,1)); %getting date vector form date number
doyvec = dayofyear(datemat(:,1)',datemat(:,2)',datemat(:,3)'); %get DOY
days_in_prev_months = [0 31 59 90 120 151 181 212 243 273 304 334 366];
lt_min = 7;
lt_max = 17;


LL = eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max;

%  % select data for the present month, lt and time limits
% LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;


time_ax = eej_data_2000_2010(LL,1);
if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;

%Period (Hour),Wavenumber  A(mA/m) B(mA/m) Std_Err_A Std_Err_B Amp(mA/m) Phase(Hours)

inv_data = ...
    [24   1 -99.64   2.55   0.84   0.30  99.67  11.90 ;
    12   2  37.38  -5.09   0.30   0.25  37.72  -0.52 ;
    8   3 -30.76  11.32   0.26   0.29  32.78  10.65 ;
    6   4  16.91  -8.44   0.33   0.30  18.90  -1.77 ;
    4   6 -13.12   1.16   0.31   0.31  13.17  11.66 ;
    24  -3  -9.31   4.97   0.31   0.31  10.56  10.13 ;
    12  -2   6.33  -3.41   0.32   0.32   7.19  -1.89 ;
    24   0   4.19  -5.76   0.32   0.32   7.12  -3.59 ;
    24  -1   6.55  -2.34   0.32   0.32   6.96  -1.31 ;
    24   5  -5.35  -1.94   0.32   0.32   5.69 -10.60 ;
    24  -2  -4.84  -2.54   0.32   0.32   5.47 -10.15 ;
    12   0  -4.64   2.65   0.32   0.32   5.34  10.02 ;
    8   5  -3.85   2.61   0.32   0.32   4.65   9.73 ;
    12   4   3.34  -3.14   0.32   0.32   4.59  -2.88 ;
    8   4  -3.17  -3.26   0.32   0.32   4.54  -8.94 ;
    12   3   2.13   3.33   0.32   0.32   3.96   3.83 ;
    12   1  -1.52   3.27   0.32   0.32   3.60   7.66 ;
    12  -1   3.29   1.23   0.32   0.32   3.51   1.36 ;
    6   3   1.51  -2.78   0.32   0.32   3.16  -4.10 ;
    24   4  -1.77   2.44   0.32   0.32   3.01   8.40 ;
    8   1   0.99  -2.30   0.32   0.32   2.51  -4.45 ;
    8  -1  -2.30   0.94   0.32   0.32   2.48  10.51 ;
    4   7  -0.85  -1.99   0.32   0.32   2.17  -7.55 ;
    6   6   1.62  -0.97   0.32   0.32   1.89  -2.07 ;
    8   7   1.16   1.35   0.32   0.32   1.78   3.28 ;
    24  -5   1.01   1.45   0.32   0.32   1.77   3.68 ;
    6   2   1.29   1.20   0.32   0.32   1.76   2.87 ;
    6   5   1.53   0.66   0.32   0.32   1.67   1.55 ;
    12  -5   1.23  -1.11   0.32   0.32   1.65  -2.80 ;
    12  -4  -0.83  -1.27   0.32   0.32   1.52  -8.22 ;
    8   2  -1.02   0.76   0.32   0.32   1.27   9.56 ;
    8  -4  -1.09   0.63   0.32   0.32   1.26   9.99 ;
    24  -6  -0.80   0.96   0.32   0.32   1.25   8.63
    ];

%% Create synthetic data and invert

imodel = 1;
idata = 12;

synth = inv_data(1,7)*cos(2*pi*((time_ax)*(24/inv_data(1,1)) + pi/(12/inv_data(1,8)) + lonnorm'*inv_data(1,2)));

for i = 2:idata,
    
    synth = synth + inv_data(i,7)*cos(2*pi*((time_ax)*(24/inv_data(i,1)) + pi/(12/inv_data(i,8)) + lonnorm'*inv_data(i,2)));
    
    
end;

% Add Gaussian Noise

%synth = synth + normrnd(0,0.1,length(time_ax),1);

% Add geomagnetic storms !

%synth(randi(length(time_ax),1,floor(length(time_ax)/20))) = normrnd(0,10,floor(length(time_ax)/20),1);




% Create a model

start_i = 1;

x_m2 = [cos(2*pi*((time_ax)*(24/inv_data(start_i,1)) + lonnorm'*inv_data(start_i,2))) ...
    sin(2*pi*((time_ax)*(24/inv_data(start_i,1)) + lonnorm'*inv_data(start_i,2)))];

for i = start_i+1:imodel,
    
    x_m2 = [x_m2 cos(2*pi*((time_ax)*(24/inv_data(i,1)) + lonnorm'*inv_data(i,2))) ...
        sin(2*pi*((time_ax)*(24/inv_data(i,1)) + lonnorm'*inv_data(i,2)))];
end;



[s,stats] = robustfit(x_m2, synth);
ls1 = x_m2\ synth;

icount = 1;

for i = 2:2:length(s)-1,
    
    fprintf('%3.0f %3.0f %4.2f %5.2f %5.2f %4.3f\n',inv_data(icount,1), inv_data(icount,2), sqrt(sum(s(i:i+1).^2)),pi/atan(s(i+1)/s(i)), ...
        sqrt(sum(ls1(i-1:i).^2)), max(stats.p(i:i+1)));
    
    icount = icount + 1;
end;

fprintf('\n     Done \n');


%% claculate the spectral aliasing

%load /Users/manojnair/projects/temp/dd.mat dd

icount = 1;

for i = 0:0.1:8,
    for j = -5:0.1:5,
        dd(icount,1) = i;
        dd(icount,2) = j;
        icount = icount + 1;
    end;
end;

synth = eej_data_2000_2010(LL,10);

for i = 1:length(dd),
    
    % create data
    
    %         synth =  100*cos(2*pi*((time_ax)*(dd(i,1)) + lonnorm'*dd(i,2)) + pi/4 );
    %
    %         % Add Gaussian Noise
    %
    %         synth = synth + normrnd(0,1,length(time_ax),1);
    
    for j = 1:length(dd),
        
        % create model
        
        
        x_m2 = [cos(2*pi*((time_ax)*(dd(j,1)) + lonnorm'*dd(j,2))) ...
            sin(2*pi*((time_ax)*(dd(j,1)) + lonnorm'*dd(j,2)))];
        
        
        % now fit the data
        
        [s,stats] = robustfit(x_m2, synth);
        
        
        spectra_matrix(i,j,:) = s;
        spectra_matrix_stat(i,j,:) = stats.p;
        
        if max(stats.p(2:3)) < 0.05 & sqrt(sum(s(2:3).^2)) > 10,
            
            fprintf('IN %3.1f %3.0f OUT %3.1f %3.0f AMP %4.2f %5.2f %4.3f\n',dd(i,1), dd(i,2), dd(j,1), dd(j,2), sqrt(sum(s(2:3).^2)),(180/pi)*atan(s(3)/s(2)), ...
                max(stats.p(2:3)));
            
        end;
        
    end;
end;



for i = 1:length(dd),
    for j = 1:length(dd),
        
        if max(spectra_matrix_stat(i,j,2:3) ) < 0.05,
            spectra_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
            spectra_phase(i,j) = (180/pi)*atan2(spectra_matrix(i,j,3),spectra_matrix(i,j,2));
        else,
            spectra_amp(i,j) = NaN;
            spectra_phase(i,j) = NaN;
        end;
    end;
end;



%% scan all the frequencies

synth = eej_data_2000_2010(LL,10);
i_count=1;



for i = 0:0.1:8,
    
    % create data
    
    %         synth =  100*cos(2*pi*((time_ax)*(dd(i,1)) + lonnorm'*dd(i,2)) + pi/4 );
    %
    %         % Add Gaussian Noise
    %
    %         synth = synth + normrnd(0,1,length(time_ax),1);
    
    j_count = 1;
    
    for j = -5:0.1:5,
        
        % create model
        
        
        x_m2 = [cos(2*pi*((time_ax)*(i) + lonnorm'*j)) ...
            sin(2*pi*((time_ax)*(i) + lonnorm'*j))];
        
        
        % now fit the data
        
        [s,stats] = robustfit(x_m2, synth);
        
        
        spectra_matrix(i_count,j_count,:) = s;
        spectra_matrix_stat(i_count,j_count,:) = stats.p;
        
        if max(stats.p(2:3)) < 0.05 & sqrt(sum(s(2:3).^2)) > 10,
            
            fprintf('%3.1f %3.0f AMP %4.2f %5.2f %4.3f\n',i,j, sqrt(sum(s(2:3).^2)),(180/pi)*atan(s(3)/s(2)), ...
                max(stats.p(2:3)));
            
            
        end;
        
        j_count = j_count + 1;
        
    end;
    
    i_count = i_count + 1;
end;

count = 1;

for i = 1:i_count-1,
    for j = 1:j_count-1,
        
        if max(spectra_matrix_stat(i,j,2:3) ) < 0.05,
            spectra_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
            spectra_phase(i,j) = (180/pi)*atan2(spectra_matrix(i,j,3),spectra_matrix(i,j,2));
            dd(count, 1) = i;
            dd(count, 2) = j;
            count = count + 1;
        else,
            spectra_amp(i,j) = NaN;
            spectra_phase(i,j) = NaN;
        end;
        
        % mean
        
        if max(spectra_matrix_stat(i,j,1) ) < 0.05,
            mean_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
        else,
            mean_amp(i,j) = NaN;
            
        end;
        
    end;
end;

%% Leakage of DE3
% Hermann suggested to look at the leakage of DE3

% 1) Period should run from 0 to 5.
% 2) DE3 Phase: -3/4 pi.
% 3) The recovered amplitudes look to me, as if they are up-scaled by a factor, say 2. I would need the unscaled amplitudes.
% 4) The labelling at the colour bar should have the same fond size as the axes labelling.
% 5) the labelling of the ordinate should be Frequency (1/day) and of the abscissa Wavenumber.

% 
frequency = 0:1:5; % day^-1
wavenum = -5:1:10;

% Synthetic data for DE3 (24 hours wave -3)
synth = cos(2*pi*((time_ax)*(1) + lonnorm'*-3)  + (pi*(-11/12)));%(1,-3)
%synth = cos(2*pi*((time_ax)*(0) + lonnorm'*0) ); % (0,0)
%synth = cos(2*pi*((time_ax)*(4) + lonnorm'*1) + pi/4); % (0,0)
synth = synth + normrnd(0,0.1,length(time_ax),1);


for i = 1:length(frequency),
    
    for j = 1:length(wavenum),
        
        %model
        x_m2 = [cos(2*pi*((time_ax)*(frequency(i)) + lonnorm'*wavenum(j))) ...
            sin(2*pi*((time_ax)*(frequency(i)) + lonnorm'*wavenum(j)))];
        
        %fit
        
        [s,stats] = robustfit(x_m2, synth);
        
        
        spectra_matrix(i,j,:) = s;
        spectra_matrix_stat(i,j,:) = stats.p;
        
        
    end;
    
end;

for i = 1:length(frequency),
    for j = 1:length(wavenum),
        
        if max(spectra_matrix_stat(i,j,2:3) ) < 0.05,
            spectra_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
            spectra_phase(i,j) = (180/pi)*atan2(spectra_matrix(i,j,3),spectra_matrix(i,j,2));


        else,
            spectra_amp(i,j) = NaN;
            spectra_phase(i,j) = NaN;
        end;
        
        % mean
        
        if max(spectra_matrix_stat(i,j,1) ) < 0.05,
            mean_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
        else,
            mean_amp(i,j) = NaN;
            
        end;
        
    end;
end;

%% plot Amplitude
fig1=figure(1);
image_data = spectra_amp;
LLL = isnan(image_data(:));
LL = isnan(image_data);

% image_data(LL) = -12;

imagesc(wavenum,frequency, image_data);
set(gca,'YTick',0:9);
set(gca,'FontSize',16);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5']);
title('Synthetic signal');
ylabel('Frequency (day^-1)');
set(gca,'XTick',wavenum);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5']);
xlabel('Wave number');
%colormap(map);
h = colorbar;set(h,'FontSize',16);   ylabel(h, 'Amplitude','FontSize',16);


[xlbl, ylbl] = meshgrid(wavenum, frequency);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','w',...
    'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',16,'FontWeight','bold');

%set(fig1, 'PaperPositionMode', 'manual');
set(fig1,'Position',[672   190   825   266]);

set(fig1, 'PaperPositionMode', 'auto');
 saveas(fig1, '/Users/manojnair/projects/tides/plots/Synthetic_test_DE3_Amp_new','png');

 %plot phase
fig2=figure(2);
image_data = spectra_phase * (24/360);
LLL = isnan(image_data(:));
LL = isnan(image_data);

image_data(LL) = -12;

imagesc(wavenum,frequency, image_data);
set(gca,'FontSize',16);
set(gca,'YTick',0:9);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5']);
title('Synthetic signal');
ylabel('Frequency (day^-1)');
set(gca,'XTick',wavenum);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5']);
xlabel('Wave number');

%colormap(map);
h = colorbar;set(h,'FontSize',16);   ylabel(h, 'Phase hours','FontSize',16);


[xlbl, ylbl] = meshgrid(wavenum, frequency);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','w',...
    'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',16,'FontWeight','bold');

%set(fig2, 'PaperPositionMode', 'manual');
set(fig2,'Position',[672   190   825   266]);
set(gcf, 'PaperPositionMode', 'auto');

 saveas(fig2, '/Users/manojnair/projects/tides/plots/Synthetic_test_DE3_Phase_new','png');


 %% Leakage of DW1
% Hermann suggested to look at the leakage of DW1

% 1) Period should run from 0 to 5.
% 2) DW1 Phase: 12 HOUR
% 3) The recovered amplitudes look to me, as if they are up-scaled by a factor, say 2. I would need the unscaled amplitudes.
% 4) The labelling at the colour bar should have the same fond size as the axes labelling.
% 5) the labelling of the ordinate should be Frequency (1/day) and of the abscissa Wavenumber.
% - The frequency scale should go from 1 to 6
% - the wavenumber from 0 to 10
% - the numbers at the colour bar should be larger
% - the heading should read ' Synthetic signal'

frequency = 0:1:6; % day^-1
wavenum = 0:1:10;

% Synthetic data for DW1 (24 hours wave 1)
synth = cos(2*pi*((time_ax)*(1) + lonnorm'*1)  + (pi*(11/12)));%(1,1)
%synth = cos(2*pi*((time_ax)*(0) + lonnorm'*0) ); % (0,0)
%synth = cos(2*pi*((time_ax)*(4) + lonnorm'*1) + pi/4); % (0,0)
synth = synth + normrnd(0,0.1,length(time_ax),1);


for i = 1:length(frequency),
    
    for j = 1:length(wavenum),
        
        %model
        x_m2 = [cos(2*pi*((time_ax)*(frequency(i)) + lonnorm'*wavenum(j))) ...
            sin(2*pi*((time_ax)*(frequency(i)) + lonnorm'*wavenum(j)))];
        
        %fit
        
        [s,stats] = robustfit(x_m2, synth);
        
        
        spectra_matrix(i,j,:) = s;
        spectra_matrix_stat(i,j,:) = stats.p;
        
        
    end;
    
end;

for i = 1:length(frequency),
    for j = 1:length(wavenum),
        
        if max(spectra_matrix_stat(i,j,2:3) ) < 0.05,
            spectra_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
            spectra_phase(i,j) = (180/pi)*atan2(spectra_matrix(i,j,3),spectra_matrix(i,j,2));


        else,
            spectra_amp(i,j) = NaN;
            spectra_phase(i,j) = NaN;
        end;
        
        % mean
        
        if max(spectra_matrix_stat(i,j,1) ) < 0.05,
            mean_amp(i,j) = sqrt(sum(spectra_matrix(i,j,2:3).^2));
        else,
            mean_amp(i,j) = NaN;
            
        end;
        
    end;
end;

%% plot Amplitude
fig1=figure(1);
image_data = spectra_amp;
LLL = isnan(image_data(:));
LL = isnan(image_data);

% image_data(LL) = -12;

imagesc(wavenum,frequency, image_data);
set(gca,'YTick',0:9);
set(gca,'FontSize',16);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6']);
title('Synthetic signal');
ylabel('Frequency (day^-1)');
set(gca,'XTick',wavenum);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6']);
xlabel('Wave number');
%colormap(map);
h = colorbar;set(h,'FontSize',16);  ylabel(h, 'Amplitude','FontSize',16);


[xlbl, ylbl] = meshgrid(wavenum, frequency);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','w',...
    'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',16,'FontWeight','bold');

%set(fig1, 'PaperPositionMode', 'manual');
set(fig1,'Position',[672   190   825   266]);

set(fig1, 'PaperPositionMode', 'auto');
 saveas(fig1, '/Users/manojnair/projects/tides/plots/Synthetic_test_DW1_Amp_new','png');

 %plot phase
fig2=figure(2);
image_data = spectra_phase * (24/360);
LLL = isnan(image_data(:));
LL = isnan(image_data);

image_data(LL) = -12;

imagesc(wavenum,frequency, image_data);
set(gca,'FontSize',16);
set(gca,'YTick',0:9);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6']);
title('Synthetic signal');
ylabel('Frequency (day^-1)');
set(gca,'XTick',wavenum);
set(gca,'YTickLabel',[' 0';' 1';' 2';' 3';' 4';' 5';' 6']);
xlabel('Wave number');

%colormap(map);
h = colorbar;set(h,'FontSize',16);  ylabel(h, 'Phase hours','FontSize',16);


[xlbl, ylbl] = meshgrid(wavenum, frequency);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','w',...
    'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',16,'FontWeight','bold');

%set(fig2, 'PaperPositionMode', 'manual');
set(fig2,'Position',[672   190   825   266]);
set(gcf, 'PaperPositionMode', 'auto');

 saveas(fig2, '/Users/manojnair/projects/tides/plots/Synthetic_test_DW1_Phase_new','png');





